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Method and Syst m for Numerically Simulating Rubb r-lik Material in 

Finite Element Analysis 

BACKGROUND OF THE INVENTION 

1 . Field of the Invention 

[0001] The present invention generally relates to method, system and 
software product used in finite element analysis, more particularly to numerical 
simulation of the structural behavior of rubber-like material under loading 
conditions. 

2. Description of the Related Art 

[0002] Finite element analysis (FEA) is a computerized method widely used in 
industry to model and solve engineering problems relating to complex systems. 
FEA derives its name from the manner in which the geometry of the object 
under consideration is specified. With the advent of the modern digital 
computer, FEA has been implemented as FEA software. Basically, the FEA 
software is provided with a model of the geometric description and the 
associated material properties at each point within the model. In this model, the 
geometry of the system under analysis is represented by solids, shells and 
beams of various sizes, which are called elements. The vertices of the elements 
are referred to as nodes. The model is comprised of a finite number of 
elements, which are assigned a material name to associate with material 
properties. The model thus represents the physical space occupied by the 
object under analysis along with its immediate surroundings. The FEA software 
then refers to a table in which the properties (e.g., stress-strain constitutive 
equation, Young's modulus, Poisson's ratio, thermo-conductivity) of each 
material type are tabulated. Additionally, the conditions at the boundary of the 
object (i.e., loadings, physical constraints, etc.) are specified. In this fashion a 
model of the object and its environment is created. 

[0003] FEA had it beginnings as a method for structural analysis, but today is 
routinely used in the design of motors, generators, magnetic resonance imaging 
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systems, aircraft engine ignition systems, circuit breakers and transformers, to 
name but a few; its techniques are used to analyze stress, temperature, 
molecular structure, electromagnetic fields, current, car crash, metal stamping, 
physical forces, etc. in all sorts of physical systems. It has become a standard 
part of the design cycle for numerous products which are not easily analyzed by 
other methods. 

[0004] FEA is becoming increasingly popular with automobile manufacturers 
for optimizing both the aerodynamic performance and structural integrity of 
vehicles. Similarly, aircraft manufacturers rely upon FEA to predict airplane 
performance long before the first prototype is built. Rational design of 
semiconductor electronic devices is possible with Finite Element Analysis of the 
electrodynamics, diffusion, and thermodynamics involved in this situation. FEA 
is utilized to characterize ocean currents and distribution of contaminants. FEA 
is being applied increasingly to analysis of the production and performance of 
such consumer goods as ovens, blenders, lighting facilities and many plastic 
products. In fact, FEA has been employed in as many diverse fields as can be 
brought to mind, including plastics mold design, modeling of nuclear reactors, 
analysis of the spot welding process, microwave antenna design, simulating of 
car crash and biomedical applications such as the design of prosthetic limbs. In 
short, FEA is utilized to expedite design, maximize productivity and efficiency, 
and optimize product performance in virtually every stratum of light and heavy 
industry. This often occurs long before the first prototype is ever developed, 
[0005] The finite element analysis method is described in detail by Thomas J. 
R. Hughes in *The Finite Element Method" (1987), published by Prentice-Hall, 
Inc., New Jersey. Generally, FEA begins by generating a finite element model 
of a system. In this model, a subject structure is reduced into a number of node 
points which are connected together to form finite elements. Governing 
equations of motion are written in a discrete form, where the motions of each 
node point are the unknown part of the solution. A simulated load or other 
influence is applied to the system and the resulting effect is analyzed using well 
known mathematical methods. 
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[0006] Two such mathematical methods are the implicit finite element analysis 
("the implicit method") and the explicit finite element analysis ("the explicit 
method"). Both methods are used to solve transient dynamic equations of 
motion and thus obtain an equilibrium solution to the equations. The methods 
march from time t through a discrete time interval or time step At. On the one 
hand, the explicit method integrates the equations of motion using the known 
quantity at time t, therefore it does not require factorization of a matrix. On the 
other hand, the implicit method solves the equilibrium conditions at time t+At. 
[0007] One of the most challenging FEA tasks is to simulate an impact event 
such as car crash. The highly non-linear behavior of the structural materials 
must be numerically simulated accurately, realistically and efficiently. The 
implicit method has been used for many years. A number of materials such as 
steel, aluminum, and rubber that used in an automobile must be included in 
such FEA software. Many components (e.g., engine and sub-frame mount 
bushings, rubber based assembly adhesives) used in an automobile are made 
of rubber, the simulation of the structural responses of njbber-like material 
becomes very important for the overall accuracy of an analysis. 
[0008] To capture the structural behavior of rubber-like material in a FEA, the 
following conditions are assumed and simulated: 1) nearly incompressible; 2) 
hyperelastic under quasi-static loading; and 3) highly rate dependent or 
viscoelastic under dynamic loading. Traditionally the hyperelastic response of 
rubber-like material has been numerically simulated using the Ogden strain 
energy function W. A brief summary of the Ogden strain energy function and 
corresponding engineering/nominal stress oo and true stress a functions are 
listed in FIG. 1. More details of the Ogden energy function is described by R.W. 
Ogden in "Non-linear Elastic Deformations" (1984), published by Ellis Honvood 
Limited, United Kingdom, which is incorporated herein by reference in its 
entirety. A number of commercially available FEA software includes these 
approaches to simulate rubber-like material. For example, LS-DYNA, a general 
purpose three-dimensional non-linear large-deformation FEA software from 
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Livermore Software Technology Corporation, is capable of simulating rubber- 
like material using the Ogden energy function. 

[0009] Today, there are a number of practical problems associated with the 
simulation of rubber-like material in FEA. Very often, engineers may only have 
axial tensile and compressive test results, which is not enough data for 
describing the required coefficients for Ogden energy function. To implement 
the Ogden energy function properly in the FEA software requires engineers to 
spend a tremendous amount of time to prepare experimental data and then 
convert them into a set of coefficients to fit a polynomial Ogden function for FEA 
software. Due to highly non-linear characteristics of this polynomial function, 
the inexactly fitted function has often resulted. This leads to a lengthy iterative 
trial-and error process of modifying the input coefficients to match the behavior 
of rubber-like material in the real world. It is therefore desirable to have a new 
method to numerically simulate rubber-like material more efficiently and 
effectively. 



SUMMARY OF THE INVENTION 

[0010] This section is for the purpose of summarizing some aspects of the 
present invention and to briefly introduce some preferred embodiments. 
Simplifications or omissions may be made to avoid obscuring the purpose of the 
section. Such simplifications or omissions are not intended to limit the scope of 
the present invention. 

[0011] The present invention discloses a system, method and software 
product to numerically simulate rubber-like material with a pragmatic approach. 
Under the assumption of material incompressibility and uniaxial loading, a 
method of calculating Ogden stress function f(A.) is developed as shown is 
FIGS. 2 and 3. The method eliminates the requirement of fitting a polynomial 
function; instead f{X) is calculated using only a set of the engineering or nominal 
stress-strain data from a uniaxial loading test for the present invention. Based 
upon the input stress-strain constitutive equation for rubber-like material, the 
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stress function f{X) is pre-computed for a range of equally spaced stretch ratios 
X. These results are stored in a computer's RAM as a lookup table and retrieved 
later during the solution phase of FEA. The advantage of the present invention 
is that no time-consuming trial-and-error process is required. The principal 
stress corresponding to a particular strain or stretch ratio can be interpolated 
rapidly from the pre-computed table directly during the lengthy solution phase. 

[0012] According to one aspect of the present invention, during the solution 
phase of FEA, the principal stresses for each integration point of the rubber 
elements is computed for each time step. In order to compute the principal 
stresses, the eigensolution of a 3x3 stretch tensor is computed first. The 
resulting eigenvalues are the three main stretch ratios for the element. The 
corresponding engineering/nominal stress is then looked up or interpolated. 
These principal stresses are converted into the true stresses and transformed 
into the global coordinate system. If the rate dependent viscoelastic effect is 
needed, the nominal stresses are looked up from a set of curves instead of a 
single curve. 

[0013] Other objects, features, and advantages of the present invention will 
become apparent upon examining the following detailed description of an 
embodiment thereof, taken in conjunction with the attached drawings. 

BRIEF DESCRIPTION OF THE DRAWINGS 

[0014] These and other features, aspects, and advantages of the present 
invention will be better understood with regard to the following description, 
appended claims, and accompanying drawings as follows: 

[0015] FIG. 1 lists the basic equations for Ogden strain energy function and 
associated stresses. The stress function f{X) is in polynomial form. 

[0016] FIG. 2 shows the derivation of the present invention based on 
incompressible material and uniaxial loading. 
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[0017] FIG. 3 shows the detailed derivation of the stress function f(X) for the 
present invention. The stress function f(X) is only dependent upon stress-strain 
curve of a rubber-like material. 

[0018] FIG. 4 shows an exemplary engineering stress-strain curve used as 
input for the present invention. 

[0019] FIG. 5 shows a process to calculate the stress function i{X) based on 
the present invention. 

[0020] FIG. 6 is an exemplary flow chart of an implementation of the present 
invention in FEA software. 

[0021] FIG. 7 depicts a block diagram of an exemplary computer, in which the 
present invention may be implemented. 
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DETAILED DESCRIPTION OF THE PREFERRED EMBODIMENTS 

[0022] In the following description, numerous specific details are set forth in 
order to provide a thorough understanding of the present invention. However, it 
will become obvious to those skilled in the art that the present invention may be 
practiced without these specific details. The descriptions and representations 
herein are the common means used by those experienced or skilled in the art to 
most effectively convey the substance of their work to others skilled in the art. 
In other instances, well-known methods, procedures, components, and circuitry 
have not been described in detail to avoid unnecessarily obscuring aspects of 
the present invention. 

[0023] Reference herein to "one embodiment" or "an embodiment" means that 
a particular feature, structure, or characteristic described in connection with the 
embodiment can be included in at least one embodiment of the invention. The 
appearances of the phrase "in one embodiment" in various places in the 
specification are not necessarily all referring to the same embodiment, nor are 
separate or alternative embodiments mutually exclusive of other embodiments. 
Further, the order of blocks in process flowcharts or diagrams representing one 
or more embodiments of the invention do not inherently indicate any particular 
order nor imply any limitations in the invention. 
[0024] To facilitate the description of the present invention, it deems 
necessary to provide definitions for some terms that will be used throughout the 
disclosure herein. It should be noted that the definitions following are to 
facilitate the understanding and describe the present invention according to an 
embodiment. The definitions may appear to include some limitations with 
respect to the embodiment, the actual meaning of the terms has applicability 
well beyond such embodiment, which can be appreciated by those skilled in the 
art: 

[0025] Loading Condition is defined as the external load acting on a structure. 
[0026] Uniaxial Load is either tension or compression. 
[0027] Longitudinal Direction is the direction, in which the uniaxial force 
applies to a material. 
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[0028] Lateral Direction is the orthogonal direction to the longitudinal direction. 

[0029] Tension is a force to pull or stretch a material. 

[0030] Compression is a force to squeeze or press a material. 

[0031] Stress o is defined as the intensity of force or force per unit area of the 

material. 

[0032] Stretch Ratio X is defined as the ratio between the deformed length and 
the original length of a material. 

[0033] Strain s is a non-dimensional quantity representing the elongation or 
shrinkage per unit length of a material under loading condition. The Strain and 
Stretch Ratio is related with the following equation: e = . 
[0034] Strain Rate is defined as the amount of strain that develops in one unit 
of time. 

[0035] Poisson's Ratio is defined as the ratio of the strain in lateral direction to 
the strain in the longitudinal direction. For rubber, Poisson's Ratio is about 0.5. 
[0036] Strain Energy is potential energy generated from the force acting on a 
material. In the elastic case, the strain energy is defined as W = 0.5*s*o- 
[0037] Deformation Gradient Matrix/Tensor is defined as the local nature of 
the deformation of a solid body. 

[0038] Principal Stretch Ratio is defined as the stretch ratio in the principal 
loading direction. 

[0039] Eigensolution is a mathematical method to determine the natural 

vibration mode shapes and amplitudes of a structure. 

[0040] Eigenvalues are the results obtained fi-om the eigensolution. 

[0041] Elastic material is a material recovers its original shape after unloading. 

[0042] Hyperelastic material is a material has one-to-one relationship between 

stress and strain for a large deformation. The stress value can be determined 

by an energy function such as Ogden stress function f{X). 

[0043] Quasistatic loading is defined as the static force dominates comparing 

to inertia and damping forces. 

[0044] Rate dependent under dynamic loading means that under dynamic 
loading the material is no longer quasistatic; for the same values of strain, 
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different values of stress may be obtained depending upon the loading velocity 
or strain rate. Rate dependent under dynamic loading and viscous material are 
interchangeably used. 

[0045] Nominal or Engineering Stress, interchangeably used herein, are force 
divided by undeformed cross section area. 

[0046] True stress is defined as force divided by current or deformed cross 
section area. 

[0047] Global coordinate system is a coordinate system used to define the 
structure in FEA. 

[0048] Embodiments of the present invention are discussed herein with 
reference to FIGS. 1-7. However, those skilled in the art will readily appreciate 
that the detailed description given herein with respect to these figures is for 
explanatory purposes as the invention extends beyond these limited 
embodiments. 

[0049] Referring now to the drawings, in which like numerals refer to like parts 
throughout several views. FIG. 1 shows the structural response functions in the 
original Ogden rubber-like material model. Ogden strain energy function W and 
true stress a\ function are listed in 110. In these functions, Xi, ?i2and 
represent the stretch ratios in three principal directions of a rubber-like material. 
J is the volumetric ratio. K is the bulk modulus, fx and a are the coefficients of 
the polynomial function f(X). These classic equations are the foundation to 
numerically simulate rubber-like material using a digital computer. Ogden 
stress function f(?i) is then defined in 120 as a polynomial function of n*^ order. 
Both engineering or nominal stress ao and true stress a are rewritten as a 
combination of f{X) in 130. It is evident that element stresses can easily be 
calculated, once 1{X) has been determined; however, the determination of the 
coefficients for the polynomial function is very difficult. 
[0050] FIG. 2 shows the derivation based on the assumption of uniaxial 
compression in longitudinal direction (?i3<1), incompressible rubber-like material 
{XA2A.3=1), and isotropic material as the symmetry of two lateral directions 
(A.i=X2). As a result of these assumptions, all three stretch ratios are rewritten in 
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the form of longitudinal direction stretch ratio X, The nominal stress go is then 
written as function of longitudinal stretch ratio X as shown in 210. In 220, X is 
replaced by eoand K is replaced by p. After equation manipulating in 220, a 
simplified nominal stress equation is shown in 230. 
[0051] The equation 230 is used as the basis for the transformation of a 
polynomial expression into a simplified deterministic expression as shown in 
FIG. 3. In 310, the expansion of the basic equation 230 is listed. It is evident 
that i{X) can be shown as a sum of an infinite series of nominal stresses. The 
formal mathematical form is listed in 320. The key significance of this equation 
320 is that there are no polynomial coefficients. Finally in 330, the nominal 
stresses and true stresses are computed with the same functions shown in 130. 
[0052] Referring now to FIG. 4, an engineering or nominal stress-strain curve 
420 is illustrated. The curve represents a set of experimental data for a rubber- 
like material under uniaxial loading condition; both tension and compression 
tests results are included. According to one embodiment of the present 
invention, the engineer needs only to prepare the experimental data such as 
this exemplary curve 420 as input to simulate rubber-like material in FEA 
software. In order to capture the stretch ratio of interest, an exemplary input 
requires a range of strain values shown in 430. The typical range of so is 
between 1 .2 and -0.8. In terms of stretch ratio X, the typical range is between 
2.2 and 0.2. If the viscous effect of rubber-like material is under consideration, 
the exemplary curve 420 may be replaced by a family of stress-strain curves 
resulting from dynamic experiments. 

[0053] With reference now to FIG. 5, because each of the next term in the 
infinite series in equation 320 depends upon the square root of the previous 
term, the order of magnitude decreases drastically from one term to the next. 
According to one embodiment of the present invention, the process 500 of 
evaluation of function f{X) is summarized in FIG. 5. At 510, f{A.) is assigned a 
value equals to X multiplied by ao(A.-1) for a given X. The stress ao value at 
strain so or X-1 is from engineering test data (e.g., the curve defined in the user 
input phase of FEA software). At 520, X is stored into a variable A^w. A new 
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variable Xnew is stored as one over square root of Xou at 530. At 540, a 
comparing test is performed for ttie absolute value of Xnew-1 being less than or 
equal to 0.01 . If the test fails, the computation of f(X) has finished, the rest of the 
terms in the infinite series is too small to affect the final result of the 
computation. If the test succeeds, another stress value at strain X,new-1 is 
multiplied by Xnewand accumulated into function f(X) at 550. At 560, the value of 
A,new is stored into Xow- The process goes back to 530 until the computation 
finishes. 

[0054] An exemplary FEA software implementation 600 of the present 
invention is shown in FIG. 6. At the input phase 610 of FEA software, at least 
one set of engineering or nominal stress-strain data for rubber material under 
uniaxial loading test is read in as a group of corresponding stress-strain pairs. If 
the loading condition is quasistatic, only one set of stress-strain curve is 
required. If the rate dependent dynamic loading condition is under 
consideration, then a family of stress-strain curves representing different strain 
rate are required. Based on the input stress-strain constitutive curve, a plurality 
of f{X) is then generated for a range of X at regular intervals at the initialization 
phase 620. These pre-computed f(A,) are stored in computer's RAM as a stress 
function table for later use in the solution phase. Next in 630, the FEA software 
calculates structural responses at each integration point of every element for 
each time step. According to one embodiment, for the element with rubber 
material, a subroutine handling rubber material is called. Within the subroutine, 
the deformation gradient matrix is gathered first at 640. For a solid element, 
either the left or the right 3x3 stretch tensors may represent the deformation 
gradient matrix. The deformation gradient matrix may be corresponded to the 
stretch ratio X or the stretch ratio squared depending upon the solution 
scheme employed. At 650, an eigensolution is solved for the deformation 
gradient tensor. The principal stretch ratio can be calculated from the resulting 
eigenvalues at 660. With these stretch ratios, at 670, the nominal stresses are 
then calculated using the corresponding f(A,) values stored in the initialization 
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phase. Because all f(?i) values of interest are pre-computed and stored as 
tabulated form, it is very computational efficient to perform a table lookup and 
interpolation. In the case of rate dependent dynamic loading, strain rate at the 
element integration point is used as basis to select two closest curves to 
determine the corresponding rate dependent stress. The true stresses are 
calculated and transformed into the global coordinate system at 690. 
[0055] With reference now to FIG. 7, a block diagram illustrates a computing 
device 700 in which the present invention may be implemented, and in which 
code or instructions implementing the processes of the present invention may 
be located. The exemplary computer system in FIG. 7 is discussed only for 
descriptive purposes. It should not be considered a limitation of the invention. 
Although the following descriptions related to a particular computer system, the 
concepts apply equally to other computer systems that are dissimilar to that 
shown in FIG. 7. 

[0056] Computer system 700 includes a processor 710 and main random 
access memory (RAM) 720 connecting to a local bus 705 through a bridge 715. 
Additional connections to local bus 705 may be made through direct component 
interconnection or through add-in boards. In the depicted example, network 
adapter 725, small computer system interface (SCSI) adapter 730, and 
expansion bus interface 735 are directly connected to local bus 705. In contrast, 
audio adapter 740, graphics adapter 745, and video adapter 750 are connected 
to local bus 705 by add-in boards inserted into expansion slots. Expansion bus 
interface 735 provides a connection for a keyboard and mouse adapter 755, 
modem 760, and additional memory 765. SCSI adapter 730 provides a 
connection for hard disk drive 770, tape drive 775, and CD-ROM drive 780. 
[0057] In order to communicate with other computer systems via a network, 
the computer system 700 connects to the network via network adapter 725. 
The network, Internet or Intranet, connects multiple network devices utilizing 
general purpose communication lines. 

[0058] Those of ordinary skill in the art will appreciate that the hardware 
shown in FIG. 7 may vary depending on the implementation. Other internal 
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hardware or peripheral devices, such as flash ROM (or equivalent nonvolatile 
memory) or optical disk drives and the like, may be used in addition to or in lieu 
of the hardware depicted in FIG, 7. Also, the processes of the present invention 
may be applied to a multiprocessor computer system. In general, Computer 
system 700 is controlled and coordinated by operating system (OS) software, 
which performs tasks such as process scheduling, memory management, 
networking and I/O services. Exemplary OS includes Linux^"^, Microsoft 
Windows™. 

[0059] Although an exemplary embodiment of invention has been disclosed, it 
will be apparent to those skilled in the art that various changes and 
modifications may be made to achieve the advantage of the invention. It will be 
obvious to those skilled in the art that some components may be substituted 
with another component providing same function. The appended claims cover 
the present invention. 
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